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1 Introduction 

Constrained Delaunay tetrahedralizations (CDTs) have many optimal prop- 
erties similar to those of Delaunay tetrahedralizations (DTs) [8, 9]. They are 
eligible structures for mesh generation. 

A major task in constructing CDTs is how to recover the domain bound- 
aries, i.e., edges and faces. At first hand, one might think it is an easy task since 
it has been addressed intensively in the literature. Unfortunately, this problem 
is far from being solved in 3D. Classical engineering algorithms [3, 12, 2] come 
with no performance guarantee. Their details (involving swaps and point in- 
sertions) are complicated to realize. On the other hand, CDT algorithms with 
theoretical guarantees are proposed [6, 7, 11]. The algorithm [11] is available 
in the program TetGen [10]. 

In this paper, we present a new algorithm for constructing CDTs for 3- 
dimensional polyhedral domains. This algorithm is based on our previous 
one [11]. A new facet recovery algorithm is proposed. It can handle facets 
which are not exactly planar - a problem ubiquitously exists in engineering 
data and the use of exact geometric predicates. We distinguish the inconsis- 
tencies between 2 and 3 dimensional constrained Delaunay simplices. These 
inconsistencies are removed by facet re- meshing and Steiner point insertions. 
This algorithm has no complex detail and is easy to implement. It can be 
extended to handle boundary consists of smoothly curved surfaces. 



2 The Algorithm 

The input is a 3D piecewise linear complex [4] (PLC) X, i.e., A' is a collection of 
polyhedra of dimensions up to 3, see Figure 1 left. We call 1- and 2-dimensional 
polyhedra of X segments and facets. The boundary complex of X is its 2- 
skeleton, which is the set of vertices, segments, and facets of X . 
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Fig. 1. Left: A 3D PLC. Right: The surface triangulation of the facets of the PLC 
with Steiner points inserted on its segments. 



Our algorithm constructs a CDT T of X. T may contain Steiner points, 
i.e., each segment and facet of X is represented by the union of a subcomplex 
of T. We call 1- and 2-dimensional simplices of these subcomplexs subsegments 
and subfaces to distinguish them from other simplices of T. 

The proposed algorithm proceeds in the increasing order of the dimensions 
of the skeletons of A'. It is divided into three phases: 

1. Create a DT Vq of the vertex set of X. 

2. Let T>\ = Vq. Recover segments of A' in a DT X>i, refine X into X' . 

3. Let T>2 =T>\. Recover facets of X' in a CDT T>2- 

The first and second phases of this algorithm are already discussed in [11 . 
Steiner points are inserted in the segments of X such that every subsegments 
are Delaunay in Pi. A symbolic perturbation technique [6] is used to ensure 
that the vertex set of V\ is in general position, i.e., no 5 points share a common 
sphere. After the second phase, the original PLC X has been refined into a 
PLC X' including Steiner points. It has been proven by Shewchuk [5] that the 
CDT of X' exists. Facet can be recovered without Steiner points. 

Shewchuk's theorem [5] is based on an assumption that all vertices of the 
facets are exactly co-planar. Unfortunately, it is commonly not the case in 
the real world data and by the use of finite precision of computer's floating 
point numbers. Therefore, a new facet recovery algorithm is developed to 
incorporate this situation. It is detailed in the following section. 



2.1 Facet recovery 

Every facet F G X' is first triangulated into a 2-dimensional CDT Tp, see 
Fig. 1 right. F is recovered in a 3-dimensional CDT T>2 when all its subfaces 
are also faces of T>2- The algorithm is summarized in Fig. 2. 

From each missing subface a one can form a missing region Q which is a 
set of subfaces such that: (i) the edges on the boundary of i? are edges of £>2, 
and (ii) the edges in the interior of Q are missing in T>2- 



A practical CDT algorithm 3 



FacetRecovery (X' , T>i) 

1. V 2 := Vi- 

2. Initialize a queue Q of all subfaces of X 

3. while Q / 0 do 



4. pop a subface a from Q; 

5. if a is missing in T>2, then 

6. form a missing region Q containing a; 

7. search a set C of crossing tetrahedra of i?; 

8. if C ^ 0 then 

9. RemeshCavity (C, i?); 

10. if J? is not recovered, then 

11. RefineRegion(A?', :D 2 , i?); 

12. endif 

13. else 

14. RemeshRegion(X> 2 , 

15. endif 

16. endif 



17. endwhile 

18. return T>2\ 

Fig. 2. The facet recovery algorithm. 



Facets with exactly co-planar vertices. In this case, the 2D CDT 

of a facet F G X' is unique viewed by all vertices of X' . If a subface a is 
missing in the current tetrahedralization T>2, then there must exist at least one 
tetrahedron (called crossing tetrahedra) in T>2 whose interior intersect with a. 
The set of crossing tetrahedra of Q forms a cavity C inside the space of X>2- The 
boundaries of C are faces in T>2. The sub-routine RemeshCavity presented 
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11] replaces the set of crossing tetrahedra by a set of new tetrahedra 
conforming to both the boundaries of C and Q. This process needs no Steiner 
points and T>2 is updated into a new CDT conforming to Q. 

Facets with not exactly co-planar vertices. In this case, there are two 
types of inconsistencies: (1) the vertices of X' may not agree on a unique 2D 
CDT of a facet F, and (2) even they agree on a unique 2D CDT, it may not 
be the current triangulation of F. Type- (2) was due to the pre-computation 
of the facet CDT assumes that the vertices of F are co-planar. To avoid these 
inconsistencies in advance would need expensive computations. The proposed 
algorithm try to handle them during the execution of the algorithm. 

A type-^ inconsistency is detected when we found no crossing tetrahe- 
dron for a missing region Q. It implies that Q can be re-meshed by a set of 
faces in T>2. The RemeshRegion sub-routine searches the set of faces in T>2 
conforming to Q and replaces the old set of subfaces in Q by the new set. As a 
result, the facet triangulation has been corrected by the constrained Delaunay 
faces of T>2. This process needs no Steiner points. 

The type-(l) inconsistency is difficult to be detected earlier. It could be 
signaled when the RemeshCavity is failed to recover Q. A typical failure 
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caused by this inconsistency is that one found that a quadrilateral abed in Q 
can not be recovered since the two triangulations of abed in the two halfspaces 
of Q are not the same. This is due to the non-coplanarity of the quadrilateral 
abed. To resolve this case, Steiner points are added into Q. We choose to 
split an edge of one of the not recovered subfaces in Q. This is done in the 
RefineRegion sub-routine. 

Note that by adding only one Steiner point may not be sufficient to achieve 
a new CDT. For an example, if this point is inserted on a subsegment of A", 
it may encroaches upon other subsegments of X' . Therefore, some tetrahedra 
containing these subsegments may not be constrained Delaunay anymore. A 
solution to solve this problem is to continue the insertion of Steiner points 
until no subsegment of X' is encroached. The RefineRegion sub-routine 
may insert several Steiner points into T>2, and update T>2 into a new CDT 
including these Steiner points. 

2.2 Termination 

The main challenge in the proof of the termination is the additions of Steiner 
points in the facet recovery phase. A key lemma we need is: The RefineRe- 
gion sub-routine returns a new CDT (with added Steiner points) of the set 
of recovered subfaces. Note that this sub-routine repairs the edges of the facet 
triangulations which are not constrained Delaunay in 3 dimensions. Once such 
edges are removed (by Steiner points), the facet triangulations can be recov- 
ered by either RemeshCavity or RemeshRegion. A complete proof of this 
lemma is our future work. 



3 Examples and Discussions 

This algorithm has been implemented in the latest version of the program 
TetGen [10] (version 1.4.3). An example is shown in Fig. 3. Comparing with 
the old algorithm, this algorithm outperforms it in both of its robustness and 
speed. Moreover, the new algorithm adds less Steiner points. One reason is 
due to the number of segments of the input is reduced when adjacent not 
exactly co-planar facets are merged into one facet. 

In the future, the termination and correctness of this algorithm need to be 
proven. The analysis of this algorithm is not complete yet. A number of ques- 
tions arise: How to incrementally update a CDT? How many Steiner points 
are needed in this process? A very interesting extension of this algorithm is 
to consider a piecewise smooth complex [1] as input. 
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Fig. 3. Example. Left: a input PLC containing non-coplanar facets. Middle: the 
surface mesh of the generated CDT. Right: a cut view of the CDT. 
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